Uniform approximation of barrier penetration in phase space 
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A method to approximate transmission probabilities for a nonseparable multidimensional barrier 
is applied to a waveguide model. The method uses complex barrier-crossing orbits to represent 
reaction probabilities in phase space and is uniform in the sense that it applies at and above a 
threshold energy at which classical reaction switches on. Above this threshold the geometry of the 
classically reacting region of phase space is clearly reflected in the quantum representation. Two 
versions of the approximation are applied. A harmonic version which uses dynamics linearised 
around an instanton orbit is valid only near threshold but is easy to use. A more accurate and more 
widely applicable version using nonlinear dynamics is also described. 



I. INTRODUCTION 

Semiclassical approaches to multidimensional tun- 
nelling lead to very interesting problems in complexi- 
fied classical dynamics, often with incompletely under- 
stood solutions. For example, recent work in references 
IH S S Q IE 111 nas shown that nontrivial geometrical 
structure such as complex homoclinic intersections have 
an important role to play in multidimensional barrier 
penetration and that even complex chaos can be relevant. 
Given the difficulty inherent in a systematic treatment 
of multidimensional tunnelling as a result of such issues, 
it is perhaps surprising that a relatively simple descrip- 
tion can be given of barrier penetration at a critical en- 
ergy where classically allowed transmission mechanisms 
turn on and where primitive semiclassical approximations 
must be replaced by somewhat more complicated uniform 
ones. 

An approach which achieves this has been proposed 
in references |7j and |8( and in this paper we apply the 
method explicitly to a model waveguide problem. The 
model is chosen to be rather simple so that fully quantum 
calculations are easy to perform accurately for purposes 
of comparison. We emphasise, however, that semiclas- 
sical aspects of the calculation are as easily applied to 
other problems, provided the topology is similar, and pro- 
vide a description, for example, of collinear atom-diatom 
reactions. For that reason we use the terminology of 
chemical reactions in this paper and equate the proba- 
bility of transmission with a probability of reaction. In 
fact the approach we describe here provides a natural 
means of visualising the quantum scattering problem in 
phase space and as such shows an interesting connection 
with classical transition state theories of chemical reac- 
tion. These classical theories have recently been of inter- 
est because the periodic-orbit dividing surface (PODS) 
construction [jj, Hfl has been generalised to arbitrary di- 
mensions using the construction of normally hyperbolic 
invariant manifolds (NHIMS) [3 13 13 M, 111 The 
classical constructions emerge naturally in our scmiclassi- 
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cal approximation and we note that even though the illus- 
trations offered here are in two degrees of freedom there 
are straightforward generalisations to higher dimensional 
problems where the full generality of the NHIM construc- 
tion comes into play. 

The approximation we use can be stated very simply 
as an abstract operator equation but for explicit illustra- 
tion we present results in phase space, using the Wigner- 
Weyl calculus. In particular, we define a Weyl symbol 
of a transmission matrix which represents, in an aver- 
aged sense, a reaction probability as a function of phase 
space. We find that above threshold the support of this 
Weyl symbol closely mimics the shape of the classically 
reacting region but the Weyl symbol itself also incorpo- 
rates tunnelling and other quantum effects. Notation and 
details for this construction are set out in section UTI 

We implement two versions of the theory. First, a 
harmonic approximation derived in is applied in sec- 
tion IIIII whose classical input consists simply of an in- 
stanton orbit, along with its action and monodromy ma- 
trix. This approximation works when the classically re- 
acting region is a small neighbourhood of the initial con- 
dition for the instanton orbit and is harmonic in the sense 
that it uses linearised dynamics generated by an ellip- 
tic quadratic Hamiltonian on the Poincare section. This 
harmonic version covers the threshold case where classi- 
cal reaction switches on as a function of energy and is 
relatively easy to apply. It fails however when the energy 
is too far above threshold and the classically reacting re- 
gion is too large to be adequately described by linearised 
dynamics. 

A semiclassical approximation that is more accurate 
and has greater range has been derived in Q and this 
is applied in section llVl This version uses fully nonlin- 
ear dynamics to extend further from the instanton orbit. 
Like the harmonic version it can be stated quite simply 
as an abstract operator equation but its practical im- 
plementation is more difficult. Difficulty arises primarily 
because we must invert an operator constructed semiclas- 
sically as an evolution operator and this inversion can- 
not at present be achieved in closed form. In this paper 
we achieve that inversion using numerical methods and 
while this aspect of the approach needs further work to 
provide an appealing semiclassical method, we can verify 
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unambiguously that the nonlinear version of the theory is 
capable of describing the quantum transmission problem 
very accurately (see Figure |HJ. 



We conclude this section by outlining how these re- 
sults relate to existing work. The basic formalism here 
of relating scattering to complex barrier-crossin g or bits 
goes back to the work of Miller and coworkers |l6l Il7j 
on the classical 5-matrix. Our intent is to describe a 
simple uniform extension of this approach which applies 
at the boundary of classical reaction where the fate of 
classical orbits changes discontinuously. What allows us 
to make progress is that we do not directly describe the 
S'-matrix but instead consider a transmission matrix de- 
rived from it which gives probabilities rather than ampli- 
tides. The advantage of this problem is that contribut- 
ing orbits at the boundary of classical reaction depend 
smoothly on initial conditions, despite the singular na- 
ture of real orbits there Q , and give simple semiclassical 
expressions. This uniformisation is similar to established 
results relating one-dimensional transmission probabili- 
ties [l8l Il9l I20I ] or the cumulative reaction probability 
|2l| to sums over multiple barrier crossings but includes 
information about how the probabilities depend on the 
incoming state. It is different however from uniform ap- 
proximations of the scattering operator such as described 
in [22J which describe explicit matrix elements. These 
are uniform with respect to variation of quantum num- 
bers whereas our approach treats the scattering operator 
abstractly and is uniform with respect to energy and in 
phase space. 
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FIG. 1: Contours of the model potential are shown for the 
case Q — 1 and A = fj. = 1/2. The continuous line shows 
a real trajectory which is an extension of the complex peri- 
odic orbit (or instanton) used in the next section to construct 
the simplest semiclassical approximation. The dashed curve 
shows the trajectory defined by the decoupled dynamics of 
Voa(y) to which it asymptotes. 



II. REPRESENTING THE SCATTERING 
MATRIX IN PHASE SPACE 

In the following sections we will develop semiclassical 
approximations for representations of the scattering ma- 
trix in phase space. In this section we illustrate these rep- 
resentations numerically using a two-dimensional waveg- 
uide, which serves as a simplified model of a collincar 
atom-diatom collision. This model consists of a particle 
of unit mass moving in the potential 

V(x, y) = sech 2 x + ^0J 2 (x)(y ~ a{x)f , (1) 



Direct approximation of the scattering problem by 
complex trajectories has recently been examined in jj, 
in the context of nonintegrable systems. It has been 
found there that, while intuitively one might expect tun- 
nelling processes to be dominated by short complex orbits 
which cross the barrier directly, the dominant complex 
orbits can have have a surprisingly nontrivial topology in 
the deep tunnelling regime. It has even been found that 
the dominant complex orbits may be chaotic in related 
treatments of quantum propagation 0, 0, 0- We also 
find evidence of the "fringed tunnelling" characteristic of 
such mechanisms in our fully quantum solutions but the 
theory we outline is intended to cover only the immedi- 
ate vicinity of the reacting region where direct tunnelling 
mechanisms are dominant. In the deep tunnelling regime 
our uniform results revert to standard primitive approxi- 
mations and we should in principle be able to marry our 
approach with that of 0,0,0 It is not obvious, however, 
that a fully uniform calculation could easily be applied 
when the contributing complex orbits are more numer- 
ous and more complicated and we do not consider that 
problem explicitly here. 



where 

uj 2 (x) = fi 2 + A sech 2 x 
a(x) = \i sech 2 x. 

A very similar potential has been used in 0, . The sim- 
plicity of this model will enable us easily to obtain ac- 
curate numerical solutions, which will be useful for later 
comparison with semiclassical approximations where ex- 
ponentially small tunnelling effects are of interest. We 
emphasise, however, that none of the theory that follows 
is dependent on this simplicity and semiclassical aspects 
of the discussion can just as easily be applied to any other 
system, as long as the Hamiltonian is an analytic func- 
tion of its arguments. The essential structural features 
we assume are that the waveguide should have a sin- 
gle bottleneck separating asymptotically decoupled chan- 
nels, which we label as reactant and product channels 
respectively, and that the energy should be sufficiently 
close to threshold that recrossings of the transition state 
do not occur. 

For future reference, it will be useful to denote the 
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asymptotically decoupled potential by the symbol 

V(x,y)~V 00 (y) = ±n 2 y 2 . 

We can therefore write the asymptotic scattering states 
for this problem analytically, as solutions of a harmonic 
oscillator. 

We write the scattering matrix in block form 

S(E) = ( \ RR tRP ) . 
v ' \tpR r PP J 

where the subscripts R and P refer to reactant and prod- 
uct channels respectively. For example, the transmission 
matrix tpp maps asymptotically incoming states on the 
reactant side to asymptotically outgoing states on the 
product side. The theory we describe is for the matrix 

71(E) = t PR t PR 

rather than for the scattering matrix itself. The ma- 
trix 71(E) has an obvious physical role determining state- 
specific reaction rates. The transmission probability for 
an incoming state labelled \ip n ) on the reactant side can 
be written as a matrix element 

Pn = (V>„|£(£)|Vn> (2) 

of 71(E). In going from the scattering matrix to 71(E) 
we lose information about phase and about the distribu- 
tion of product states but, as described in 0, uniform 
semiclassical approximations for 71(E) are expected to 
be considerably simpler than those for S(E). 

Formally, we can think of 71(E) as an operator acting 
on the Hilbert space H^(E) of asymptotically propagat- 
ing states in the incoming reactant channel — in this 
context we refer to it as the reaction operator in the fol- 
lowing. We consider systems for which there are a finite 
number M of such states so Ti 1 ^ (E) is finite-dimensional 
and 71(E) can be represented by an M x M matrix. The 
Wigner-Weyl correspondence offers an alternative repre- 
sentation as a function on phase space and the theory we 
outline in the coming sections is stated in those terms. In 
the remainder of this section we describe how the connec- 
tion is made formally between the matrix representation 
and the representation in Wigner-Weyl correspondence, 
or the Weyl symbol of 71(E). In preparation for that dis- 
cussion, let us first describe the phase space on which 
these representations are defined. 

As described in p| (see also the natural classical 
analog of the Hilbert space H^(E) is a Poincare section 
S 1 ^ (E) defined by fixing the energy E and an asymptot- 
ically large value of the reaction coordinate x. We use 
(y,p y ) as canonical coordinates on T,^(E) and we may 
alternatively denote points in T,^(E) using the vector 
notation 




We define the allowed region of the y-p y plane by the 
condition p^/2 + V oa (y) < E and note that, as usual in 
this sort of correspondence, the dimension M oiH^(E) is 
approximated by the Liouville area of this region divided 
by 2nh. 

To calculate the Weyl symbol W^(C,E) of 71(E), we 
denote its individual matrix elements by R nm (E) = 
(ijj n \7i(E)\ip m ) and write 

W^iCE) = 2nhJ2 Rn m (E)W nm (C) (3) 

nm 

where W nm (C) are the Weyl symbols of the projectors 
|V'n)(V'm|/27r?i (the factors of 2ttH are to keep the notation 
consistent with standard practice for Wigner functions in 
the case n = m). For the asymptotically harmonic poten- 
tial inffl the functions W n m(C) are known analytically 
(see [24| for example). 

The Weyl symbol W^(C,E) of 71(E) provides a re- 
markably transparent means of visualising the quantum 
transmission problem and of relating the scattering ma- 
trix to the geometry of classical phase space. To illus- 
trate, we show examples of W^(C, E) in Figure[3]for the 
model potential in Q with energies at and above thresh- 
old. These results have been obtained by first comput- 
ing the scattering matrix numerically using symplectic 
integration combined with the log derivative method as 
described in [2^ and then using ©. In each case 
we see that W^(C,E) is effectively supported in a re- 
gion of Ej!j(.E), which we can identify as the quantum- 
mechanically reacting region. Indeed, by rewriting the 
reaction probability in (using standard properties of 
the Wigner-Weyl correspondence) in the form 

Pn = J W^(C,E)W nn (C)dC, (4) 

where d£ = dydp y , and interpreting the Wigner func- 
tion W nn (C) as a phase space pseudodensity, it is natural 
to identify W^(C,E) as a probabability of reaction as a 
function of phase space, albeit in an averaged sense. Al- 
though the uncertainty principle prevents us from defin- 
ing a point- wise transmission probability in phase space, 
we can construct linear combinations of incoming states 
with a fixed total energy (as in J^J below) whose Wigner 
function is supported within an area of 0(h) in H^(E) 
and the appropriate modification of @ then gives the 
reaction probability as an average of W^(C, E) over that 
support. 

For energies at or below threshold, transmission is con- 
trolled by tunnelling and W^(C,E) is supported in a 
phase space region of area 0(H), centred around an initial 
condition that leads to an optimal tunnelling route. An 
explicit semiclassical expression for VV^C, E) in this case 
will be given later and one can see for the threshold case 
in Figure |2t a) that W^(C, E) is indeed peaked around 
a single point in S^(£J). As the energy increases above 
threshold, a classically reacting region appears, initially 
centred on the orbit associated with optimal tunnelling. 
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FIG. 2: The Weyl symbol W^(£,E) is shown for the potential in Q. Cases (a) to (d) have = 1 and A = n = 1/2 and 
energies (a) E = 1.0, (b) E = 1.10, (c) E = 1.15 and (d) E = 1.6. Cases (e) and (f) have Q = 1, A = -1/2 and /i = 1 and 
energies E = 1.15 and E — 1.6 respectively. Also shown are, the boundary of the allowed region as a dashed curve and, in 
cases (b) to (f) where the energy is above threshold, the boundary of the classically reacting region. In each case a box of area 
h is shown on the bottom left. 



The boundaries of the classically reacting regions are in- 
dicated in Figures|2b) to (f ) by continuous closed curves 
(as the energy falls to the threshold case E = 1, the clas- 
sically reacting region shrinks to a point corresponding to 
the optimal tunnelling route and around which W 7 j(C, E) 
is concentrated). One can see in each case that the re- 
acting region closely matches the support of W^(C,E). 

Before describing how VVU(£, E) is approximated semi- 
classically, we should outline how the classically reacting 
regions in Figure |21 are defined. The boundary of the 
classically reacting region in full phase space is the stable 
manifold on the reacting side of a PODS in two dimen- 
sions or more generally a NHIM if higher-dimensional 
problems are treated. Extended into the incoming reac- 
tant channel, this stable manifold defines a tube, the in- 
terior of which consists of classically reacting trajectories 
and whose annular exterior in the classically allowed re- 
gion consists of trajectories which eventually return along 



the outgoing reactant channel. A representation of this 
reacting region in a Poincare section T, 1 ^ is obtained sim- 
ply by taking the a section of the tube of reacting tra- 
jectories at fixed energy and a fixed, asymptotically large 
value of the reaction coordinate x. Since the tube contin- 
ues to evolve asymptotically (according to the dynamics 
of a decoupled potential Voo(y)), the shape of a reacting 
region defined in this way will depend on the value cho- 
sen for the coordinate x. In models where the dynamics 
of the reacting region is nonlinear — or potentially even 
chaotic in problems of higher dimension — the shape of 
the reacting region will not have a limit and becomes ever 
more complicated as x is brought to infinity. To obtain 
a fixed asymptotic limit we therefore renormalise the dy- 
namics by using the decoupled evolution of the limiting 
potential Voo(y) to map the asymptotic section back to 
one corresponding to a fixed finite value of x. In making 
semiclassical comparisons the value of x used to define 
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this final section is dictated by the conventions used for 
the scattering matrix. 

In the present case the asymptotic states in terms of 
which the scattering matrix is defined are of the form, 




ipn(y), 



(5) 



where tp n (y) are the eigenfunctions of the decoupled po- 
tential Voo(y). The phases of these scattering states are 
zeroed at x = and asymptotic incoming states can 
be constructed by starting with the transverse modes at 
x = and propagating them backwards into the asymp- 
totic region of the incoming channel using dynamics de- 
fined by V oa (y). To make a comparison with the classi- 
cal picture, the renormalisation of the classical dynam- 
ics should therefore take an asymptotic Poincare section 
back to one defined by x = 0. This is the convention used 
in FigureElto compare W 7 j(C, E) with the classically re- 
acting region. Alternative phase conventions would lead 
to a reaction operator 71(E) obtained by conjugation of 
the one we define by a unitary matrix which is diagonal 
in the basis \ipn)- This conjugation makes no difference 
to the diagonal matrix elements (ifj n \7V\^ n ) but is im- 
portant for the appearance of the Weyl symbol in Eg. 
Choosing values other than x = for the reference sec- 
tion Eg would, for example, lead to a deformation of the 
Weyl symbol by the asymptotically decoupled dynamics. 

Note that this renormalisation procedure is simply 
means of interpreting a term in the phase function in 
the classical S'-matrix 01 that fixes its asymptotic value. 
By accounting for this term using a conjugation of the 
asymptotic dynamics by the mapping in decoupled dy- 
namics back to i = 0, we can incorporate everything 
about the classical S'-matrix into a single Poincare map- 
ping and present results in a more compact form, as de- 
scribed more fully in the coming sections. 



III. HARMONIC APPROXIMATION 

In this section we describe a semiclassical approxima- 
tion for 7Z(E) based on dynamics linearised around an 
optimal tunnelling orbit. Although less accurate and 
valid over a smaller range of energies than the fully non- 
linear theory described in the next section, this har- 
monic approximation captures the essential qualitative 
behaviour of 71(E) and works well in the especially inter- 
esting range of energies around threshold where classical 
reaction switches on. It is also considerably simpler to 
apply and can be expressed in closed form using easily 
obtained classical data. Note that the term "harmonic" 
here refers simply to the fact that linearised dynamics 
about a tunnelling orbit are used and has nothing to do 
with the harmonic asymptotic behaviour of the model we 
use for numerical illustration. 



A. Operator version 

A full description of the harmonic approximation to 
71(E) has been given in 7] and we refer there for de- 
tails and a derivation of the approach. Here we simply 
summarise the important points. The operator 71(E) is 
approximated by a formula 



71(E) 



?(E) 
1 + T(E) 



(6) 



where T(E) is a "tunnelling operator" constructed from 
classical data. The elements needed to compute T(E) 
are as follows. 

• A complex periodic orbit 75 (t), the "instanton", is 
found which encircles the transition state in imag- 
inary time. This orbit can be found for energies 
above and below threshold and its dynamical char- 
acteristics depend smoothly on energy there. 

• At the end of the previous section it was described 
how the boundary of the reacting region can be first 
extended arbitrarily far into the asymptotic region 
and then renormalised by mapping back to a sec- 
tion Eg at x = using decoupled dynamics, so that 
the shape remains fixed as dynamics are extended 
into the asymptotic region. An analogous renor- 
malisation, illustrated in Figure [3 is applied to 
"fE (t) so that an asymptotically fixed initial cond- 
tion for it is defined in Eg. Above threshold this 
initial condition is near the centre of the reacting 
region. 



• The imaginary action of 7_e(i) is denoted iKo = 
§ p • dq. We also denote 9 = K /h. Note that 
for energies above threshold we have Kq < while 
Kq > below threshold. 

• Linearised dynamics around 7_e(£) are characterised 
by a complex monodromy matrix W, which is rou- 
tinely determined as part of a numerical search for 
the orbit je (t) ■ The eigenvalues of W come in real 
reciprocal pairs (A,A _1 ), which we order so that 
A > 1. 

• The matrix W can be generated by using an elliptic 
quadratic Hamiltonian /i(C) for an imaginary time 
—Itq. Note that this Hamiltonian generates renor- 
malised dynamics in the section Eg and is therefore 
not simply a truncation of the full Hamiltonian in 
the transition state region. 

• Canonical coordinates (Q, P) are defined on Eg so 
that 

h(C) = f (0 2 + p 2 ) 

and we have A = e aT ° . 
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FIG. 3: The harmonic approximation to VV 7 j(C, E) is shown for energies (a) E = 1, (b) E = 1.10 and (c) E — 1.15, corresponding 
to the top row in Figure [5] 



The quantum analog of the classical generating Hamilto- 
nian is denoted by h. 

We can now write the tunnelling operator T(E) in the 
form 



(7) 



which, except for the prefactor e 6 , is an imaginary-time 
evolution operator generated by h. Since h is harmonic 
we can explicitly construct its eigenstates \<fk), with A; = 
0, 1, • • • and the corresponding eigcnsolutions of T are 

T\ifk) = Tk\(fk) 

where the eigenvalues 

T k = e- 9 A-( fe +3) (8) 

are deduced simply by exponentiating the eigenvalues of 
h. 



B. Weyl symbol 

It is shown in 7] how closed form approximations can 
be deduced for phase space representations of 71(E) as a 
result of substituting this exponentiated form for T(E) 
in JJJJ and resumming the geometric series 1Z = T — T 2 + 
T 3 — ■ ■ ■ . This leads to an integral representation 



1 r P -p(e+T h/h) 
K(E) = T . / : dp 



(9) 



for 71(E) , in which the contour C ascends just to the right 
of the imaginary axis. Standard asymptotic approaches 
to this integral, such as the method of steepest descent, 
allow explicit asymptotic approximations to be written 



for TZ in various representations, including for the Weyl 
symbol. 

These expressions are especially useful to understand 
the detailed structure of 71(E) in phase space, but for the 
purposes of computing 7Z(E) for the parameter regimes 
we consider here, it suffices to use an an eigenexpansion 



7i(E) w yy fc |<ifc)(ifci 



where 



Tk 



Tk 



(10) 



(11) 



1+Tfe 

The Weyl symbol of 7Z(E) can, for example, be written 



as 



W^(C,E) « 2nhY,rk(E)Wkk(C) 



where Wfefc(C) are the Wigner functions of the states \ifk) 
(and given analytically in |2 1| for example). The tilde 
distinguishes these Wigner functions from those of the 
basis states \tp n ) of the scattering operator, which are 
different. 

The canonical coordinates (Q, P) are centred on the 
initial condition for je (t) in and are such that for 
energies just above threshold, the classically reacting re- 
gion is circular in the (Q, P) plane. In the original coor- 
dinate system (y,p y ) these Wigner functions are trans- 
lated, squeezed and rotated so that their level curves are 
aligned with the approximately elliptical reacting region. 
The resulting approximation for W^(C,E) therefore de- 
scribes an elliptically-shaped representation of the true 
quantum transmission problem 

The harmonically approximated Weyl symbol 
W^(C,E) is illustrated in Figure [3] for three energies at 
and just above threshold in the model potential with 
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2, 




FIG. 4: The harmonic approximation is compared with exact 
results for W^(^, E) at energies at and just above threshold. 
These are obtained by sampling W^(C, E) along a horizontal 
line through the centre of the reacting region in E 1 ^ and plot- 
ting the result as a function of y. Parameters in the potential 
are as in Figures|2Ia)-(d) and|3] 



1.5r 1.5 




-2-1 1 2-2-1 1 2 



(a) (b) 



1.5 r 1.5 




-2-1 1 2-2-1 1 2 



(c) (d) 

FIG. 5: The harmonic and exact Weyl symbols are compared 
for different potentials with Q. = 1 and (a) A = — |, fx = 0; 
(b) A = 0, n = 0; (c) A = §, A* = 0; (d) A = §, fJ, = |. In each 
case the energy is chosen so that the classically reactive flux 
has the fixed value N C \(E) = 2. 



fl = 1 and A = 1/2 = /x. For comparison, illustrations of 
corresponding exact calculations can be found in the top 
row of Figure |21 In general we find that the harmonic 
approximation is in good quantitative agreement with 
exact results at and below threshold. The threshold case 
in Figure|3Ja), for example, is indistinguishable from the 
corresponding exact result in Figure EI a) at the level 
of graphical resolution used. As energy increases, the 
agreement deteriorates so that noticable differences are 
visible when E — 1.15 (harmonic approximation in Fig- 
ure |3{c) and exact calculation in Figure Hfc)). It should 
be emphasised, however, that even then, the harmonic 
approximation captures the essential qualitative features 
ofW^(CS). 

In order to make a closer comparison between exact 
and harmonic results, we show one-dimensional sections 
through the Weyl symbol in Figure In each case 
W^(C,E) is sampled along a horizontal line through 
the centre of the reacting region in and plotted as 
a function of the y coordinate. There is good quan- 
titative agreement in cases (a) and (b) where the en- 
rgy is at and just above threshold. At higher energies 
the harmonic approximation captures the support of the 
quantum-mechanically reacting region well but details of 
the Wigner function do not match at the centre of the 
reacting region. It should be remarked, however, that 
oscillations in the Weyl symbol are sensitive to nonlocal 
changes in phase space and discrepencies at the centre 
of the reacting region may not have a strong effect on 
averaged reaction probabilities as expressed in 

Similar one-dimensional sections are shown in Figure^] 
which illustrate the harmonic approximation for different 
parameter sets. In each of these the energy is chosen so 



that the cumulative reactive flux 

Ncl( - E) = 2^h f PODS p ' dq = nV^Tx 

is fixed (at the value 2). Figures GJa) , (b) and (c) show 
cases where /i = and the potential has a symmetry in y. 
Figure [SJ a) has a negative value of A = —1/2 for which 
an adiabatic approximation assuming fast transverse dy- 
namics in the barrier region would not be expected to 
apply. Figure E^b) is the separable case A = and in 
Figure OJc) we have A = 1/2. Note that separability 
does not confer a particular computational advantage in 
this approach, nor does it lead to particularly better ac- 
curacy of the approximation. In Figure^d), an example 
is shown in which ji = 1/2 and the potential is neither 
separable nor symmetric in y. The approximation works 
less well in that case. This is not unexpected because 
corrections to the harmonic approximation will be quar- 
tic rather than cubic in a symmetric problem but we note 
that there is still good agreement. 

IV. NON-LINEAR APPROXIMATION 

Although the harmonic approximation captures the es- 
sential qualitative features of quantum transmission and 
works well quantitatively near threshold, we can achieve 
greater range of applicability and significantly improved 
numerical agreement if we use fully nonlinear dynamics 
around the orbit (t). The price to be paid for this 
improvement is that the resulting calculation is signifi- 
cantly more involved. The greatest impediment is that, 
although the tunnelling operator defined by nonlinear 
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evolution can be routinely approximated semiclassically, 
we do not know at present how to write semiclassical ap- 
proximations for the operator (l+T) -1 directly in terms 
of classical orbits. In this paper we simply use numerical 
inversion of the matrix representation of 1 + T . Before 
describing this procedure, it is helpful to describe how 
nonlinear calculation is incorporated in the operator T. 
This is done in section llV Al below, followed by a descrip- 
tion of the uniform calculation in section IIV Bl 



A. Primitive approximation 

At energies below threshold the imaginary action of 
the orbit 7e(£) is positive, that is 9 > 0, and the ex- 
ponential prefactor e~ e in 10 makes T small. We may 
therefore approximate the reaction operator directly by 
the tunnelling operator, giving 

which we refer to as the primitive approximation. The 
primitive approximation is easily extended beyond the 
immediate neighbourhood of 7_e(*)- Instead of letting T 
be the evolution operator corresponding to the classically 
linear evolution defined by W, as we did in the previous 
section, we let it be the quantum version of a nonlinear 
map in Eg. 

Initial conditions near 7e(£) in Eg can be followed 
over a sequence of time evolutions similar to those of 
je (t) itself until they return to Eg, defining a surface- 
of-section mapping which we denote by 

J- . Ij r — * Zj r . 

As with conventional return maps, T defines a canonical 
transformation on Eg, except that it is complex, in gen- 
eral taking real initial conditions to complex images. De- 
spite this complexity, the evolution has a quantum ana- 
log as an evolution operator, and this is the tunnelling 
operator T . 

With suitable modifications to take account of the 
complexity of the mapping, standard semiclassical ap- 
proximations that are applied to evolution operators, 
such as the Van Vleck formula, can be used to approxi- 
mate T. Here we focus on an approximation derived in 
[2?| for the Weyl symbol of an operator (see also [28|h 
For the Weyl symbol of the operator T we write 



-A(C,E)/H 



y/det {W AB + I)/2 ' 



(12) 



where A and Wab are calculated from a midpoint orbit 
C A — > £ B which is defined by the conditions 



c 

Cb 



1 



(CA + Ci 




-4 



x 



FIG. 6: The projection onto real configuration space of a 
typical midpoint orbit is shown. This orbit can be obtained 
by continuously deforming the initial conditions for the orbit 
shown in Figure The dotted segment corresponds to evo- 
lution over the imaginary part of the time contour and the 
dashed segment corresponds to renormalisation by the un- 
coupled dynamics. The symmetry £ s = C,* A means that the 
real part of the orbit shown here is self-retracing, with the 
orbit turning back on itself at the end of the dotted segment. 



That is, the point £ at which the Weyl symbol is to be 
evaluated is the midpoint of C, A an( i Cb > where £ A evolves 
into Cb under the return map. The exponent A(C, E) is 
such that 



iA(C,E) 



C B 



p-dq-p y (y B - Da) 



and the matrix Wab is a linearisation the map T around 
the orbit Ca ^ Cb- 

It can be shown ||| that the complex conjugate of the 
map T is its inverse, T* = J 7 ^ 1 , and from this a num- 
ber of important symmetries follow which guarantee that 
VUf (Cj E) is a real- valued function on Eg that is peaked 
around the initial condition for 7_e(£), which we denote by 
£ in the following. On a formal level, the real-valuedness 
of W.f (C, E) follows from the observation that T is Her- 
mitian which, as discussed in 28] , is a quantum analog of 
the property T* = T~ x . It is instructive, however, to see 
how the real-valuedness of the semiclassical approxima- 
tion to W-f (C, E) follows directly from the symmetries of 
the midpoint orbit. 

First we note that, given a midpoint orbit £ A — > £ B for 
a real- valued £, then £* B — > £* A is also a midpoint orbit 
(for the same £). This can be seen by conjugating the 
relations in JTSJ and using C B = [{^(Ca)}* = F*(C*a) = 
T~ (Ca) ^0 deduce that £* A — JF(Cg). It turns out in 
fact that these two midpoint orbits coincide, so 



and 



Cb = Ca 



ReC B =C = ReC 



(13) 



This is easily confirmed for the linearised map (replac- 
ing T by multiplication by W and using W* = IF -1 ) 
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FIG. 7: For an energy E = 0.985 just below threshold, a 
comparison is given on a logarithmic scale between the ex- 
act Weyl symbol (continuous curve), the nonlinear primitive 
approximation of Equation il 1 2li (heavy dashed curve) and 
the harmonic primitive approximation of Equation 1 1 H (light 
dashed curve). The parameters here are f2 = 1, A = —1/2, 
and n = 1. The nonlinear primitive result works well over a 
large region of phase space but does not capture the oscilla- 
tions in the tail of the exact calculation, which may well be 
a signature of nonintegrable complex-dynamical effects of the 
type described in 0, 

and therefore holds for the nonlinear map if C, is close 
enough to £ . The condition £ B = £ A can only be 
violated if a bifurcation is encountered and a more de- 
tailed analysis shows that this corresponds to the condi- 
tion det (Wab + I) = 0, which would lead to a caustic 
in (|12fl . We will assume in this paper that no such caus- 
tics are encountered in the region of S 1 ^ which dominates 
reaction. 

An example of a full trajectory corresponding to a typ- 
ical midpoint orbit is illustrated in Figure^ Because the 
initial conditions are complex, the coordinates of the tra- 
jectory are generically complex over its length, even along 
the segments which have been obtained by deformation of 
the real segments of 7e(£)- A consequence of the symme- 
try Cb = Caj however, is that the time contours can be 
chosen so that the second half of the trajectory reverses 
the complex conjugate of the first half. A projection onto 
real configuration space, for example, is self retracing. 
We find as a result that the action is purely imaginary 
and the exponent A(C, E) is a positive real number. We 
also find that W\ B — W A b an d because det Was = 1 
this means that the amplitude term det (Wab+I) in 
is real (and positive). Therefore Wf(C,E) is a positive 
real-valued function with a maximum at £ (for which 
we have A(£ , E) = K (E)). 

The harmonic approximation can be recovered by ex- 
panding the exponent to second order about £ and ap- 
proximating Wab by W, giving 

-e-[(2//3)tanh/3/2]r h(C)/fi 

ssm ■ (14) 

A comparison is given in Figure [7] between this approx- 



imation, the fully nonlinear approximation of l|12|) and 
exact results. Although the harmonic approximation 
works well near the maximum of the Weyl symbol, the 
fully nonlinear result works better over a larger range. 
We note however that there is qualitative deviation even 
from the nonlinear approximation in the deep tunnelling 
regime where the exact calculation shows significant os- 
cillations not captured by the harmonic or nonlinear ap- 
proximations. Similar oscillatory structure, or "fringed 
tunnelling" in the scattering matrix has been explained 
in 0- IE M on the basis of nonintegrable complex dynam- 
ics and has been shown to involve mechanisms that also 
show up in chaotic tunnelling. It seems likely that the 
oscillations in Figure \7\ have a similar origin but we have 
not preformed a detailed analysis. It will be an inter- 
esting problem in the future to combine the inherently 
nonintegrable mechanism in 0, 0, [|| with the uniform 
approximations illustrated here. 

B. Uniform approximation 

Although we now have an explicit closed-form semi- 
classical approximation for 7", no equivalent result is cur- 
rently available for the uniformisation 7/(1 +T) because 
inversion of the operator 1 + 7" cannot be done simply. 
In this paper we simply adopt a hybrid approach which 
combines semiclassical approximation of 7" with numeri- 
cal inversion of 1 + 7". Although not a fully semiclassical 
method, this will allow us to verify that © gives an ac- 
curate reproduction of quantum transmission. 

We first represent T as a matrix in the same asymp- 
totic basis \ip n ) as used for the scattering matrix. We 
denote individual matrix elements by T nm = (ip n \'I'\i^m) 
and compute them using the Wigner-Weyl calculus by 
writing 

T nm {E)= fw nm (C)W f (C,E)dC (15) 

and approximating Wf(C,E) using (|12fl . We emphasise 
that while T is almost diagonal in the basis \<Pk) of eigen- 
states of the generating Hamiltonian h, the same is not 
true in the basis \ip n ) unless the potential is separable. 
The integral is performed numerically and the resulting 
M x M matrix for 7"/ (1+7"), whose elements are denoted 
R^ m (E), is also computed numerically. This integration 
is not difficult since a grid on E 1 ^ is easily filled by using 
Newton integration to step the midpoint orbit, starting 
with the known solution corresponding to 7e(i) at £ 
(for which Ca = Co = Cb)- Once the elements R s £ m {E) 
are known, the Weyl symbol for 7~/(l +7") is obtained 
by replacing R nm (E) with R s £ m (E) in 

We find excellent agreement between the semiclassi- 
cally computed Weyl symbol W 7 ^(C, E) and the exact re- 
sult. A nonlinear version of Figure is indistinguishable 
from the exact results shown in the top row of Figure [21 
and therefore not shown. Instead we compare in Figure^] 
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FIG. 8: The uniform nonlinear approximation is compared 
with the exact Weyl symbol for a series of energies and the 
same parameters as used in Figure 0] It is difficult the distin- 
guish the exact and approximate results at the level of graph- 
ical resolution used, but in each case the difference multiplied 
by 10 is shown underneath as a dotted curve. 



horizontal slices of the Weyl symbol through the reacting 
region, in the same manner as in Figure The potential 
used is the same as in Figure 0] and the energies treated 
extend somewhat higher above the threshold. The exact 
and semiclassical results cannot be distinguished at the 
resolution used, so the difference scaled by a factor of 
10 is also shown. We find similarly good agreement for 
other parameter sets we ahve investigated and we note 
that the quality of the approximation does not require 
special features of the classical dynamics such as sym- 
metry, separability or adiabatic separation of transverse 
from reaction degrees of freedom. 

We should remark that the current hybrid implemen- 
tation of the nonlinear calculation is cumbersome and is 
harder to apply further above the barrier where the larger 
region of integration demands that we extend the mid- 
point trajectory deeper into complex phase space. We 
have not, for example, reproduced the results on the sec- 
ond row of Fig. [fusing this method. The purpose of this 
calculation is to show that the nonlinear uniform result 
derived in || provides an accurate description of TZ(E) 
in the model considered and that the method therefore 
deserves further exploration. Even though numerical in- 
version was used in applying the formalism, it is built en- 
tirely on a semiclassical approximation for the tunnelling 
operator T(E) and we expect that any subsequent fully 
semiclassical implementation will be equally accurate. 

We also remark that the current hybrid method is the- 
oretically clumsy and obscures somewhat the deeper con- 
nections between the quantum results and the underlying 
classical geometry. For example, it would be especially 
interesting to characterise the behaviour of W^(C, E) at 
the boundary of the classically reacting region where tra- 



jectories approach the PODS (or NHIM in higher dimen- 
sions) along its stable manifold and where the classical 
reaction probability drops sharply from 1 to 0. Such an 
analytical approximation was found in for the har- 
monic version in which W 7 ^(C, E) is approximated as an 
integral of the Airy function near the boundary of clas- 
sical reaction. Investigation is currently underway into 
a method to derive similar results in the nonlinear case 
on the basis of generating T as in equation J7J), but with 
an anharmonic generator h computed using classical nor- 
mal form theory. Ultimately it should be possible to de- 
scribe explicitly how the quantum reaction probability 
varies across the boundary of the classically reacting re- 
gion in terms of trajectories which approach the complex- 
ified PODS along its stable manifold and evolve along it 
before returning to the asymptotic Poincare section. 



V. CONCLUSIONS 

We have successfully treated quantum transmission 
across a multidimensional barrier using uniform semi- 
classical approximation. The method applies generically 
around a threshold energy and does not rely on specific 
features of the classical dynamics such as separability or 
the existence of action angle variables. In its fully nonlin- 
ear incarnation the method gives an accurate description 
of reaction probabilities in phase space and makes a strik- 
ing connection between the quantum scattering problem 
and the geometry of classical reaction. We expect that 
it will work equally well in higher-dimensional problems, 
even in cases where the incoming states are chaotic in 
the transverse dynamics. 

Although we have shown that the fully nonlinear ver- 
sion works well, in doing so we have resorted to numerical 
methods which are not in the spirit of semiclassical ap- 
proximation. The theory therefore needs further develop- 
ment in order to achieve a fully semiclassical description 
of the emerging reacting region. One promising approach 
which is currently under investigation is to use classical 
normal form theory to generate dynamics around the or- 
bit je (t) using a nonlinear extension of the generator 
h(C). Many of the explicit analytical approaches used in 
the harmonic case might then be adapted to the nonlinear 
approximation. In particular this is expected to produce 
a detailed analytical description of the Weyl symbol at 
the boundary of the reacting region which calls on in- 
trinsic geometrical features of (the stable manifold of) 
the NHIM. 

A second aspect of the calculation which deserves fur- 
ther attention is the treatment of rotational degrees of 
freedom in fully three-dimensional models of chemical re- 
action. Although at one level this is simply a question 
of applying the results here individually to symmetry- 
reduced phase spaces for given angular momentum quan- 
tum numbers, there are interesting and nontrivial prob- 
lems in describing the quantum-classical correspondence 
compactly in operator form. This is an especially in- 
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